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I. INTRODUCTION 

The development of plasma fluid models that incorpo- 
rate finite-Larmor-radius (FLR) corrections has a rich 
history in plasma physics (see Ref. [l| and references 
therein and Ref. [2j]). Plasma fluid models offer the great 
advantage of simple tractability of a few fields with ex- 
cellent small-scale resolution of space-time scales. One 
of the disadvantages of fluid models over kinetic models, 
however, is that specific features of the orbital particle 
dynamics are lost by averaging over particle momentum 
(or velocity) space. 

Various fluid models have, thus, been built by incorpo- 
rating one or more of the following kinetic effects: stan- 
dard FLR effects associated with the transformation from 
particle phase-space dynamics to guiding-center phase- 
space dynamics [!, 0, [a, SLQlj wave-particle resonances 
(e.g., Landau damping) [8|, Q; and separate contribu- 
tions from magnetically-trapped and untrapped particles 
in tokamak geometry [l(| lllj . The Hamiltonian struc- 
ture of reduced magnetohydrodynamic equations [l2| has 
also led to the development of sophisticated reduced fluid 
models [3, [l3[ that retain a Hamiltonian structure as well 
as preserve important conservation laws. In particular, it 
is often stressed in the derivation of reduced fluid mod- 
els that the energy-conservation property (important for 
accurate numerical implementation) represents an impor- 
tant property to be preserved through reduction. The ex- 
istence of a Lagrangian variational principle from which 
reduced nonlinear fluid equations are derived guarantees 
the derivation of an exact energy conservation law by the 



Noether method. 

The purpose of the present work is to (i) generalize 



the previous work by Brizard [14j . (ii) provide a new 
interpretation for the reduced-fluid equations presented 
in Ref. [Tij in terms of nonlinear FLR corrections, and 
(iii) demonstrate the properties of these equations with a 
linear dispersion relation (to uncover the types of waves 
they describe) as well as linear and nonlinear numerical 
simulations in straight and magnetic dipole geometries 
to investigate energy conservation. 



A. Perturbed electric and magnetic fields 



First, the perturbed electric and magnetic fields 



b CM|| 

E = - V$ * 

c dt 



(1) 
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used in the present work are expressed in terms of the 
perturbed scalar potential $ and the perturbed vector 
potential Am bo, where we assume that the perturbed 
magnetic field is perpendicular to the (time- independent) 
background magnetic field B = B n b n . 

Note that the formula for the perturbed magnetic field 
used in Ref. [141 ] was not in general divergenceless (un- 
less V X bo = 0, i.e., bo is the gradient of a scalar), 
but Eq. is explicitly divergenceless. Because of this 
modification, some of the equations in this paper differ 
from those of Ref. In particular, while magnetic- 

curvature effects were absent from the parallel Ampere 
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equation in Ref. [lj], these effects are taken into ac- 
count in the present work [l5| . Note also that the condi- 
tion b • = implies the assumption b • V X bo = 
(bo/So) * V X Bo = 0, i.e., we assume that no back- 
ground current Jug cx bo • V X Bo flows along the back- 
ground field lines (e.g., magnetic dipole geometry). This 
assumption is adopted here in order to simplify our dis- 
cussion of nonlinear FLR effects and we note that the 
condition bo • V X Bo ^ can easily be restored in our 
nonlinear reduced fluid equations (see Ref. [l6| for exam- 
pie). 



B. Standard and nonlinear FLR corrections 

Second, the new interpretation for the reduced-fluid 
equations of Ref. [l4[ is given in terms of nonlinear FLR 
corrections that arise from nonlinear gyrokinetic the- 
ory [It]]. Standard FLR corrections are associated with 
the transformation from particle phase space to guiding- 
center phase space. In particular, the gyroradius vector 



T gc x x 



(3) 



is interpreted as the displacement between the represen- 
tation of the particle position in guiding-center phase 
space Tg^x and the guiding-center position x (here, Tj^. 1 
denotes the push-forward operator associated with the 
guiding-center phase-space transformation; see Appendix 
[A] for further details). 

Standard FLR corrections in plasma reduced-fluid for- 
malism are based on the fact that a guiding-center par- 
ticle feels electromagnetic fields that are averaged over a 
gyration period. Hence, a guiding-center particle (with 
mass m and charge q) feels the gyroangle-averaged po- 
tential 



(*(x + P0)> = $(x) 



pBp 
2mOo 



V x $(x) 



, (4) 



where, by definition, the gyroangle-average of the lowest- 
order gyroradius p = p H vanishes (i.e., (p ) = 0) 

and the lowest-order nonvanishing FLR correction in- 
volves the gyroangle-averaged dyadic product {p p } = 
h Po I J. j where f2 = qB / mc denotes the (signed) gy- 

rofrequency, po — [2p,Bo/ '(toSI^)] 1 / 2 , and Ij_ = I — bobo. 
The ordering of the standard FLR correction in Eq. (TJJ 
involves the factor A = k\p\, where k± = \k±\ de- 
notes the perpendicular wavenumber. In Fourier k-space 
(where V_l — > ikj_)j the expansion Q may be summed 
up as (<&) — > Jo(A)pk, where Jo(A) is the zeroth-order 
Bessel function (ESQ- 

In nonlinear gyrokinetic theory [r?1 |. the asymptotic 
elimination of the gyroangle dependence reintroduced 
by the perturbation of guiding-center dynamics by low- 
frequency electromagnetic fluctuations is carried by the 
transformation from guiding-center phase space to gy- 
rocenter phase space. The combination of the guiding- 



center and gyrocenter phase-space transformations intro- 
duces the gyrocenter displacement 



t: 



1 gc < 



(5) 



where T gy 1 (T gc 1 x) denotes the representation of the 
particle position in gyrocenter phase space and x de- 
notes the gyrocenter position (here, Tgy 1 denotes the 
push-forward operator associated with the gyrocenter 
phase-space transformation) . The difference between the 
guiding-center position x and the gyrocenter position x 
is directly proportional to the perturbed electromagnetic 
fields d)-© [see Eq. © below]. The lowest-order term 
in the gyrocenter displacement vector p gy = p 1 + ■ ■ ■ is 
expressed in terms of the gyrocenter phase-space coordi- 
nates (5c, Vu , p, C) as [TtJ 



Pi 



{Si, x + p }gc 
g f dSi dp 
mc \ <9£ dp 

bo dSi cb 
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qB 



dSi dp 
dp d( 



(6) 



Here, { , } gc denotes the guiding-center Poisson bracket 
(higher-order terms involving spatial gradients of the gy- 
roradius Po have been omitted) , and the gyrocenter scalar 
field is given by the first-order FLR expression 



Si = 



n d( 



E 



c 



b X B± 



(7) 



where Ex = — Vj_$ and higher-order (standard) FLR 
terms have been omitted. The gyroangle-average of the 
gyrocenter displacement vector (jSJ) yields 



{Pi) 



({Si, Pols 



_q_ d_ I dSi 

mc dp \° <9£ 



BqQq 



bo X B J 



(8) 



where the guiding-center contribution (p ) and the gyro- 
center contribution (Si) combine to yield a nonvanishing 
gyroangle-averaged gyrocenter displacement. 

The addition of the gyrocenter displacement vector (J6j) 
to the standard gyroradius vector p thus means that the 
gyrocenter particle now experiences an averaged poten- 
tial 



($(x + p + Pl )) = $(x) 



• V$(x) 



pB =2 

V ± $(x) 



2mf2 2 



> (9) 



that combines the standard (guiding-center) and nonlin- 
ear (gyrocenter) FLR corrections. Note that the rela- 
tive importance of nonlinear FLR effects relative to the 
standard FLR effects is represented by the ratio q$/T 
and, thus, nonlinear FLR effects dominate in the cold- 
ion-fluid limit when q$/T > 1 [TJ, [TJ], i.e., when the 
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linear E x B velocity (normalized to the thermal speed 
Wth = y/T/m) is large enough so that it satisfies the con- 
dition |uE|/t>th 3* fcj_Pth- In order to focus our attention 
on the new nonlinear FLR corrections that appear in 
the reduced fluid equations of Ref. [14J, however, we set 
k±pth = in the present work and postpone our discus- 
sion of standard FLR corrections in reduced fluid models 
to future work (see Ref. [l6| for standard FLR corrections 
to a reduced electrostatic fluid model). 

The combination of the guiding-center and gyro- 
center phase-space transformations gives the relation 
Tgy 1 (Tg^x) = x + p + p 1 between the particle position 
x and the gyrocenter position x. Since the gyroangle- 
average of the gyrocenter-particle displacement 

(t-^t-V) - s) = < Pl ) jt o (io) 

does not vanish, it leads to the well-known polarization 
and magnetization effects in the gyrokinetic Maxwell's 
equations pjj (see Appendix [5] for details concerning 
the dynamical reduction of the Vlasov equation by Lie- 
transform method [2p]|). 

By averaging the gyrocenter displacement vector (|8|) 
over the gyrocenter Vlasov distribution in gyrocenter mo- 
mentum space, we obtain the reduced-fluid displacement 
vector 

where um denotes the gyrocenter-fluid parallel velocity, 
tiii H±/Bo represents the magnetic flutter velocity, and 
the linear perturbed E x B velocity is 

u E ee |^XV ± $ = E x xf^. (12) 

Using Eq. (jTTJ) , we define the effective potentials 

$ p = $ - p ± -E ± ] 

\ > (13) 
A\\ p ee A\\ - b -p ± X B ± J 

which include nonlinear finite-Larmor-radius (NFLR) 
corrections that are quite distinct from the standard FLR 
corrections; see Appendix [B] for further details concern- 
ing a physical interpretation of the reduced-fluid displace- 
ment (fTTj) as well as a Lie-transform derivation of the ef- 
fective potentials (fl~3"|) . The purpose of the present work 
is to investigate how polarization and magnetization ef- 
fects manifest themselves in reduced fluid equations self- 
consistently derived from a Lagrangian variational prin- 
ciple. 

The reduced gyrocenter-fluid moments (n, u,pj_,p\\) 
used in the present work, obtained as moments of the 
reduced Vlasov equation 0, are also expressed 



in terms of the physical (phys) fluid moments and the 
reduced-fluid displacement (fTTj) . According to Appendix 
|A1 the reduced-fluid density and parallel velocity U|| = 

bo • u are 

"-phys = n - V - (npjj H 1 

\ > (14) 

where higher-order nonlinear FLR effects arc ignored. 
The reduced-fluid perpendicular and parallel pressures 
are expressed in terms of similar nonlinear FLR expan- 
sions. Lastly, we note that the treatment of higher-order 
gyrocenter-fluid moments (e.g., heat fluxes) is presently 
outside of the sco pe o f a variational formulation (see Ap- 
pendix A of Ref. [14j for additional comments). 

C. Energy conservation properties 

Third, the energy conservation properties of our non- 
linear reduced fluid equations are guaranteed by the use 
of a variational principle. By using the Noether method 
[2l| . the local energy conservation law 

l+V.S-0, (15) 

is derived with explicit expressions for the energy den- 
sity £ and the energy-density flux S. To verify energy 
conservation in our numerical simulations, we decompose 
the energy density £ = J2i &t m terms of the compo- 
nents £i (with corresponding energy-flux decomposition 
S ee Si) and track the time evolution of each volume- 
integrated component JSj ee J v £i d 3 x in terms of the 
energy-transfer equations: 

d -§ = E Q«> (i 6 ) 

3 

where Qij = — Qji denotes the (antisymmetric) energy 
transfer between components i and j (such that total en- 
ergy conservation is guaranteed by ^ Qij = 0) and 
boundary conditions are chosen such that the surface 
terms § gv Si • r\dA vanish in Eq. (|16p . 

D. Organization 

The remainder of the paper is organized as follows. 
In Sec. |TtJ we review the derivation of the reduced fluid 
equations and the self-consistent Maxwell's equations of 
Ref. [3] by the Lagrangian variational method. The 
derivation differs from the one presented in Ref. [l4[ in 
our treatment of the perturbation magnetic field @ . In 
Sec. IIII1 we derive the exact energy conservation law (fTS"|) 
by Noether's method. In Sec. IIV1 we rearrange the re- 
duced parallel- force equation derived in Sec.|TTJto display 
the nonlinear FLR effects explicitly based on Eq. (|13p . 
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In Sec. El we present a linear dispersion relation for a 
homogeneous two-fluid isotropic plasma. In Sec. IVIi we 
present linear and nonlinear numerical simulations us- 
ing our reduced-fluid equations in straight and magnetic 
dipole geometries and we summarize our work in Sec. lVIll 
In Appendix [XJ we present a summary of the general 
foundations of polarization and magnetization effects as- 
sociated with dynamical reduction in plasma physics pre- 
sented in Ref. [20j and, in Appendix [B] we present a sim- 
ple interpretation of the reduced-fluid displacement (fTT]) 
as well as a simple derivation of the effective potentials 

(USD. 



II. REDUCED-FLUID EQUATIONS 

The nonlinear finite-beta reduced-fluid equations de- 
rived by Brizard [l4| were obtained from a variational 
principle 



0. 



(17) 



with the perturbed nonlinear E x B velocity 



cbo 

R2 
n 



(E x X Bx) 



Bx 

B ' 



(22) 



The parallel reduced- fluid velocity U\\ represents motion 
along the perturbed magnetic field lines, the effective po- 
tential <&* includes the zero-Larmor-radius (ZLR) gyroki- 
netic electrostatic correction, and the effective potential 
At includes the parallel nonlinear E x B velocity 



At 



bo 



Bx 

B 



■ An b 



We will see below that these definitions simplify the equa- 
tions of motion obtained from the variational principle 

(HU). 

It is useful to introduce the following nonlinear finite- 
Larmor-radius (FLR) identities involving the reduced- 
fluid displacement vector (fTTj) and the definitions (fT3")) : 



where the Lagrangian density (sum over species is im- 
plied unless otherwise noted) is 



-f|Ej 



IBI 



qn ■ 



An 



ran 
~2~ 







.Bx\ 


u \\ 


(b -\ 






B J 



A|,boJ (18) 



(qn$ + V). 



Here, V = p± + ^ p\\ = \ Tr(P) is the trace of the Chew- 
Goldberger-Low (CGL) pressure tensor 



P = P||b b + px(I-b bo). 



(19) 



Note that the Lagrangian density (JTSJ) contains the linear 
E x B velocity u E and the linear magnetic-flutter velocity 
[22l [23 | u\\ R±/Bo explicitly, while all other terms have 
standard interpretations. In addition, we note that the 



replacement |E| 



E, 



in the first term in Eq. (fT5)) 



removes the parallel displacement current (dtEu) in the 
parallel Ampere equation and the total magnetic field 
is B = Bo + Bx, where the background magnetic field 
Bo = V X Ao is assumed to be time-independent. 

Following Ref. [l4| , the Lagrangian density (fT5)) is also 
written as 



1 

8^ 



E 



IBI 



+ — U \\ - (<i n * 

where the effective fields are 

b n 



qn 
c 

V) 



(a -u + u\\ A\ \ 

(20) 



Bx/Sc 



= U|| b 



<7$* = q§ - m \u E \ 2 /2 
{q/c)A\ = (q/c)A\\ - mV\\ 



(21) 



q$* + -U( 



m ., 



Kn, 



•At 



mu|| b 2 



A 



\p 



where the second-order Hamiltonian term 



k p = ^ng| Pj 



(23) 
(24) 

(25) 



is interpreted as a low-frequency ponderomotive Hamil- 
tonian term. Indeed, in nonlinear gyrokinetic theory [l7| . 
the second-order gyrocenter Hamiltonian term 



\ ({Si, {S 1} Ho} }) 



O 2 
"0 



2 B dfj, 




m 

T 



yields, upon averaging with respect to the gyrocenter 
Vlasov distribution over gyrocenter momentum space, 
the relation 



— n 



IPX 



2 = K 



(26) 



where we have neglected thermal effects. The form of this 
ponderomotive Hamiltonian is similar to the magnetic- 
moment Hamiltonian fiB = imfi 2 , (|Pol 2 ) (i- e -i the 
guiding-center "ponderomotive" Hamiltonian). It is also 
similar to the high-frequency ponderomotive potential 
[24| K w = i mu}' 2 WpIjW, where || • ■ ■ || denotes eikonal- 
phase averaging and the high-frequency eikonal displace- 
ment p w = — (e/moj' 2 ) (E + v/c X B) is expressed in 
terms of the high-frequency wave electric and magnetic 
fields (with u>' = u> — k- v). 
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A. Reduced polarization and magnetization vectors 



From the Lagrangian density (| 18[> . we define the follow- 
ing reduced-fluid polarization and magnetization vectors 



Mi = 



SL 


_ Ex 


SE ± 


47T 


dL 




dB ± 


An 



= £«npj_, (27) 

J_> v— v / till — X 

— = "£qn p x X-ibo) ! (28) 



which are expressed in terms of the reduced-fluid dis- 
placement vector (fTTj) . We note that the perpendicular 
magnetization vector (|28[) is expressed in terms of the 
moving-dipole contribution only [25j. Note also that the 
relation 

Y, nK P = J( E ^- p i + Bi-MJ (29) 



exemplifies the K~x Theorem [26J, |27|, |28j , which leads to 
the expression 



IBI 



where we have defined the macroscopic electromagnetic 
fields 



D = E ± + 47rPj 
H = B-4ttMi 



(30) 



The reduced Maxwell's equation are, thus, expressed as 
V-D = 4tt 9™> ( 31 ) 

and 



f SD 



V X H — = 4tt > -nu, 



(32) 



where the right sides represent the total gyrocenter 
charge density and the total gyrocenter current density, 
respectively. 

We may thus rewrite the Lagrangian (|20|) as 



L = 



sfO 



E, r - iBi 



U till 

V + qn(A -~ + -±A llp -$ p! . 



" I y u \\ - K < 



(33) 



where nonlinear FLR corrections JT3J) and the reduced- 
fluid ponderomotive (|26p are shown explicitly. 

In the absence of nonlinear FLR effects (i.e., p ± =0), 
the Lagrangian density ([33]) reverts back to the standard 
Lagrangian density for a guiding-center plasma fluid. 



B. Dynamical constraints 

The variational principle (fT7j) does not treat the fields 
(n, u,pj_,p||; E, B) as independent variational fields. In- 
stead, the Eulerian variations (8n, Su, Sp±, Spn) are ex- 
pressed in terms of the virtual fluid displacement £ while 



the Eulerian variations (<5E, SB) are expressed in terms 
of the potential variations (£)<£>, SA) subject to constraint 
equations. 

The constraint equations for the Eulerian fluid- 
moment variations (8n, 8u, 5p±, 8p\\) are the continuity 
equation for each reduced-fluid species 



nK p = -^(Ej_-D - B-H), 



Sn 
St 



= - V-(nu), 



(34) 



and the CGL pressure equations 
dp± 



dt 
S^| 

at 



V-(p±u) - p± (I-bobo): Vu, (35) 
V-(p||u) - 2p,| bobo : Vu, (36) 



where the higher-order heat-flux moments are omitted 
here (but were considered in Ref. [Hj]). 

The Eulerian fluid-moment variations Sr) a = 
{Sn, Sp±, Sp\i) are defined in terms of the relations 

^ s 2h>(^r At )> (37) 

where the virtual fluid displacement £ is defined as 



lim (uAt) 

At^o 



According to Eqs. l(5g] l -Po ]) and Eq. (07]), the Eulerian 
variations 8rj a are 



8n 
8V 




(38) 



and the Eulerian variation 8u = Au — £ • Vu of the fluid 
velocity is defined in terms of the Lagrangian variation 

d£ S£ 

Note that the Eulerian variations 5n and <5u satisfy the 
constraint St<5n + V • (Sn u + n 8u) = 0, as expected. 

The constraint equations for the Eulerian variations 
(£E, SB) are 

1 SB 

VXE+- — = 0andV-B = 0. 

c St 

The Eulerian variations for the electromagnetic fields E 
and B: 

SE = - V<5$ — and <SB = V X SA (39) 

c at 

are expressed in terms of the variations and 5 A. 

We express the variation of the Lagrangian (12T))) in 
terms of (£, (5$, SA) through the relations (f3"8"]) and 
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SL = -£ 

+ V-P 



d fdL\ „ / dL 



Vu- 7T 
au 



/a£ 

5A -{8A 



dL 
dn a 
I d_dC 
cdtffE 



, , dC „ 9£ 

5$ 1- V 

9$ <9E 



V X 



ac 

9B 



IF 



v-sr, 



where the space-time-divergence components 
dL 



6A 
ST 



du 

( dL 
\ du 



- <^A • — - , 

c 9E 5 



<9n a 



„ 9L dL 

w ot + 5Ax aB 



(40) 



(41) 



(42) 



do not play a role in the least-action principle (|17jl but, 
instead, play a crucial role in the derivation of exact con- 
servation laws (in the next Section). 



C. Reduced-fluid equation of motion 

The stationarity of the action (fl7j) with respect to £ 
yields the Euler-Poincare equation for the reduced-fluid 
velocity u: 



d fdL\ _ / dL 

= at Us +v ' u aS 



dL\ 
V-P - n V— 

an / 



— 

<9u 



(43) 



which becomes the reduced-fluid equation of motion 
dul 



= ran b 



where 



II 
dt 



qn[F,* 



X B* + V-P*, (44) 



6 2 and V • P* = V • P + \ mnWUh and 



| L< cfcJJLVJ. v - I V "I I "2 f/Wfc V 

we introduced the effective electric field 



E* 



= _ V $* - — 1 
c at 

= E + -(~V|u,| 2 



and the effective magnetic field 



(45) 



B 



V x 



A 



B + V x 



mc 
mc 

q 



Mil 6 2 ) b 



U|| b 2 ) b 



(46) 



We clearly see that the definitions (u|,$*,^4|) used to 
write Eq. (|44p in simple form actually hide all the non- 
linear corrections. 



The reduced-fluid equation (|44|) may be decomposed 
into two separate equations. First, the cross-product of 
Eq. ([44]) with bo yields the reduced-fluid velocity 



u = um b* 
where 



cb 
qn B 



- X (V-P* + <jnV$*), (47) 



BJ = b -B* and b* = (48) 



Note that, under the assumption bo • V X bo = 0, we find 
B* = B and b* = B*/B . Next, the dot-product of 
Eq. (|4~4"|) with b* yields the reduced-fluid parallel equation 
of motion 



du*, 

mn—± = b* • (qnE* -V-P* 
dt w 



(49) 



Once again all nonlinear corrections are hidden in the 
definitions of the effective fields (it*, <&*, which may 
prevent us from arriving at a clear interpretation of these 
nonlinear terms. The simplicity of Eq. (|49[) . however, 
points toward some underlying principle behind the orga- 
nization of the nonlinear terms. We postpone providing 
our new interpretation of Eq. (I49p in terms of nonlinear 
FLR corrections until Sec. ITVl 



D. Reduced Maxwell's equations 

We now use our variational principle (|17[) . based on the 
variation (|40"|) . to derive the reduced Maxwell's equations, 
which exhibit the important polarization and magnetiza- 
tion effects that provided the motivation for nonlinear 
gyrokinetic theory [l7l |. 

The stationarity of the action J SL d 3 x = with re- 
spect to (5<I> yields the Euler-Poincare equation 



ac „ dc 

a¥ +V 'aE = °' 
which becomes the reduced Poisson equation 
V-E, 



(50) 



4?r 



^] ?[n - V-(npx)] = Pphys- (51) 



Here, the physical charge density p p hy S is expressed as 
the sum of the reduced charge density (VJ en) and the 
polarization density (— V • P_l). 

The stationarity of the action (|17jl with respect to 5 A 
yields the Euler-Poincare equation 



0, 



dC 1 d_ dC dC 
dA + cdtdF, + X dB 
which becomes the reduced Maxwell equation 



(52) 



B 



^_ 1 > v — > 11 1 a / e 1 

V X — = y -nu + - — I— - + P 

47r ^ c c dt \ 4n 



V x Mi 



I phys • 



(53) 
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Here, the physical charge current J p h ys is expressed as the 
sum of the reduced charge current (^~J enu), the polar- 
ization current (dP±/dt), and the magnetization current 
(cV X Mi). The parallel component of Eq. ([53]) yields 
the reduced parallel- Ampere equation 



b n - V X 



Bi 

47T 



Mi 



E 



qn 



(54) 



Note that in a straight and uniform background magnetic 
field (V X bo = 0), the reduced parallel- Ampere equation 
(15*11) becomes 



V 2 ± A 



4tt 



E 



qn 

c 



V-i 



'-II Pi 



which appeared in Ref. [141 ]. The present work, however, 
makes use of the more general equation (I53|) . In the spe- 
cial case of zero equilibrium current (such as for a dipole 
field), b • V X B/(4tt) = - V • [V ± (A 11 /B )BI}/(AttB ). 

III. ENERGY CONSERVATION LAWS 

One of the great advantages of using a variational prin- 
ciple to derive exact or reduced dynamical equations re- 
sides in the fact that these self-consistent equations are 
guaranteed to possess important conservation laws (e.g., 
energy-momentum or wave-action) . This is especially im- 
portant for reduced fluid equations that are derived by 
imposing an approximation scheme based on space-time- 
scale orderings on exact fluid equations. 

When the three Euler-Poincare equations P"3"j) . ("5"0"]) . 
and lf5"2")) are taken into account, the variation of the La- 
grangian (|40| reduces to the Noether equation 



SL 



dSA 
~~dt 



V-ST. 



(55) 



The Noether equation (|55| may be used to derive con- 
servation laws for the reduced fluid equations (|M|) - (f!*!t5)) . 
(JUJ), (J5TJ>, and with the reduced-fluid velocity u 

given by Eq. ([44]) . 



A. Local energy conservation law 



We use the Noether equation (|55|) to derive a local 
energy conservation law associated with time-translation 
symmetry (i — > t + St), where 



£ = 
5$ = 
SA = 
SL = 



uSt 

Std t § 
5td t A 
Stdth 



= cSt (E + V$) 



(56) 



Upon rearran ging terms and performing some gauge can- 
cellations [201 |21|, we obtain the local energy conserva- 
tion law (115[) . where the energy density is 



£ 



dh 
du 



dh 



E- 



dh 



(57) 



and the energy-density flux is 
ulu-g) + IP 



drf 



c E X 



dh 
dB 



dL_ 
OA 



(58) 



By substituting derivatives of the reduced-fluid La- 
grangian (|20|). Eqs. ([57]) and ([58]) become 



£ 



8tt 



(|Ej 



IBI 



E 



D 

47T 



_ mn / 9 

V + — ( U \\ - \"e\ 



(|Ej 



IBI 



and 



mn 



S = - — E_l X H 

47T 

mn 



B_ 

(P 



U 



VI) u 

2^ 



(59) 



(60) 



which are identical to those presented by Brizard [T3| if 
we take into account that the perturbed magnetic field 
Bj_ is now divergenceless. 

B. Global energy conservation law 

By combining the CGL pressure equations (|35[) and 
()36|) . we obtain the evolution equation for the internal 
(pressure) energy density 

^ = -V- (uP + P-u) + u-(V-P), (61) 
ot 

where the energy-flux term appears in Eq. (|60p . By 
combining the continuity equation (|34[) with the reduced 
parallel-acceleration equation |49|) . we obtain the evo- 
lution equation for the parallel kinetic energy density 
(which includes motion along perturbed magnetic field 
lines) 



d /mn 
dt VT~ 



/mn n\ 



u-(V-P) 



db 



mi - | qE* - mC7||b —J , (62) 

where the energy-flux term appears in Eq. (|60[) . Lastly, 
by dotting the reduced Maxwell equation ("53|) with Ej_, 
we obtain the evolution equation for the electromagnetic 
energy density 

d_ 

di 



E 



D 



1 



4?r 



8 71 



(|Ej 



IBI 



— nu • 



— Ej_ X H - u 

47T 

qE - mU\\b — 



(63) 
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where the energy- flux term appears in Eq. (|60p . 

By adding the evolution equations ([6"Tj) - (j6"3")) , all trans- 
fer terms (i.e., terms that are not exact divergences) can- 
cel each other exactly and we recover the local energy 
conservation law (|15[) . where the energy density £ and 
energy-density flux S are given by Eqs. ([59]) and ([60]) . re- 
spectively. We note that, if we label the evolution equa- 
tions ([5T )) -([55 j) as dEi/dt + V • S, = J2j^i Hi with * = 1 
(internal energy), i — 2 (parallel kinetic energy), and 
« = 3 (electromagnetic energy), then the antisymmet- 
ric energy-transfer density matrix q^ = — qji has the 
nonzero components 



912 = u-(V-P), 
q 23 = nu- [q~E* 



m f/iib, 



db 
dt 



(64) 
(65) 







1) 









The energy-transfer equations (|16p therefore become 



Ql2 

Qu + Q23 ) , (66) 

-Q23 



where Qij = J v qij d 3 x = — Qji denotes a component of 
the volume-integrated energy-transfer matrix. Equation 
(f6"6")) shows the importance of the effective parallel kinetic 
energy (E%) in the transfer processes with the internal 
(pressure) energy (E\) and the electromagnetic energy 
(E 3 ). 



NFLR EFFECTS IN REDUCED-FLUID 
PARALLEL DYNAMICS 



IV. 



The reduced-fluid parallel equation of motion (|49|) 
gives the time evolution of the effective parallel velocity 
field ut = it|| b 2 = U|| (1 + |B^| 2 /_Bq) in terms of the effec- 
tive fields (|2ip. Equation (|4"9")l is written in a form where 
the convective part u-Vun is hiding on the right side. 
This equation may be written in a more standard form 
that brings out the NFLR corrections (fT3"|) explicitly. In 
order to facilitate our interpretation of this equation, it 
is preferable to write it in a way that explicitly displays 
the total time derivative, du\\/dt = dfUu +u • Vitii, of the 
parallel reduced- fluid velocity un. 

First, using the nonlinear FLR identities (|2"3")l and 
Eq. (|49]) becomes 



n ^- (mu\\ + - Ay^j = b* ■ (F p — ran U||Vit||) , 

(67) 

where we introduced the NFLR-corrected force density 



F p = - V • P - n V (q $p + Kp) 
= - V • P p - qnV$ p . 

The effective CGL-pressure force density 



V-P 



V-P 



n VK n 



(68) 



(69) 



now contains the low-frequency ponderomotive force den- 
sity nVKp. This ponderomotive correction appears in 
complete analogy with the high-frequency ponderomo- 
tive force density [2!| that appears in reduced fluid mod- 
els US El- 
Next, we write the magnetic vector (|48|) as 



Vtt|| X 



where 



b* = b + -!L V X b 

"0 



B 



Bo 



(70) 



(71) 



includes the standard guiding-center curvature term 
(u||/f2o) V X bo and the NFLR-corrected perturbed mag- 
netic field 



Bi P = V X (An, b 



UP-..;- ( 72 ) 
so that the reduced-fluid velocity (|4T|) may be written as 



uu b* 



^P x b„ 



(73) 



9 ' mn f^o 

After rearranging terms in Eq. (149|) and using the iden- 
tity 



-mnunb* • Vu|| = 



-mn || b* • Vit|| 



bo 



= — mnU'Viiii — X Vun • F p 



-mn u • Vuii 



(b!-b*).F„ 



we finally obtain 

du\\ 



dt p V c " dt 

= b* p ■ (qn E p - V • Pp) , 

where the NFLR-corrected electric field is 

b <9A||p 



E„ 



(74) 



(75) 



c dt 

It is immediately clear that in the absence of nonlinear 
FLR and ponderomotive effects (i.e., p ± = 0) and omit- 
ting background magnetic curvature, the reduced-fluid 
parallel equation of motion ([74]) reverts back to the stan- 
dard parallel equation 

"'" ' b„ - r± ] . (gnE - V-P), 



dt 



B_ 

So 



where the total time derivative d/dt — d/dt + u- V is 
expressed in terms of the guiding-center fluid velocity 



u = u 



B_ 

B 



cb 
qnB Q 



X (gnV$ + V • P) , 



which includes the standard E x B and magnetic-flutter 
convective nonlinearities. With the nonlinear FLR and 
ponderomotive effects retained, Eq. (|74p expresses the 
parallel momentum equation in the displaced (by pj_) 
frame of the gyrocenters. 
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TABLE I: Parameter definitions and values for the linear dispersion relation (Fig.[T}, the linear simulation in straight geometry 
(Fig. EJl, an d the linear and nonlinear simulations (Figs. [4] & [5} in dipole geometry. 



Symbol 



Definition 



Figure 1 



Figure 2 



Figures 4 & 5 



V 

K 

A 

fitot 



(k±c/io pl ) 2 
Va/c 
m e /mi 
3/3 s /2 
Pl+P'e 



0.01 
0.001 
0.05 
0.1 



1 

0.05 
0.01 
0.15 
0.3 



0.001 
0.02 
0.001 
0.045 
0.09 



V. LINEAR DISPERSION RELATION 

For a two-component plasma fluid in a homogeneous 
magnetic field (for which the equations of this paper are 
identical to those of Brizard [lj]), the linear dispersion 
relation of our equations is expressed (in terms of the 
symbols defined in Table |J) as 

( v -TT^:) [l ' (1+f "' ) ^« l 

+ [t vl V'' -V(% + «„/j;) + #/SJ] K = 0, (76) 



where Va = B$/ y/Aitmirii is the Alfven speed, u P i = 
\J ATTTiiqf j mi is the ion plasma frequency, (3 S = 8irp s / B 2 
is the plasma beta for species s (= i or e for ions or 
electrons), and 711 (= 3) is the ratio of specific heats for 
the parallel (ID) motion. 



The first line of Eq. (|76|) (neglecting the terms propor- 
tional to re, and with e m and e c small) yields Alfven wave 
(V ~ 1) and sound wave (V ~ f3' tot ) solutions. These 
basic waves are modified by finite re, e m , and e c . Further 
details concerning the interpretation of Eq. (jT6[) are given 
in Ref. [H, HI] . 

Figure [T] shows the normalized squared phase speed 
(defined in Table U by V) and the ratio of relative 
magnetic fluctuations to relative density fluctuations 
6B y /5n e = (5By/Bo)/(5n e /Tio) versus the normalized 
squared perpendicular wavevector re (defined in Table 
HI with k± = k x ) using the linear dispersion relation 
Eq. ([75)1 . Here, Sn e = n e — n is the perturbed gyro- 
center density rather than the perturbed physical density 
(Sriphys)e = n e — V • (n e p ±e )~n [Eq. {H}]. However, for 
the electrons, V • (n e p ±e ) is small (because p ±e cx m e ) 
so 5n e is approximately equal to the perturbed physical 
density. (On the other hand, rii is not a good approx- 
imation for the physical ion density unless the density 
perturbations are dominated by fluctuations caused by 
parallel motion, as for low beta sound waves.) Note in 
Fig. [T] that at low re, there are Alfven (solid curve) and 
sound wave (dashed curve) solutions. The Alfven solu- 
tion has V = 1 (w/fcii = Va) an d small density fluc- 
tuations relative to magnetic fluctuations. The sound 
wave solution has V = /3' tot [(ui/k\\ ) 2 equal to the squared 
sound speed 7||(Tj + T e )/mj] and small magnetic fluctua- 
tions compared to density fluctuations. Unlike the MHD 
equations, our equations have coupled parallel fluid mo- 



tion and parallel current [Eq. (|53"|) ]. so there are coupled 
magnetic and density perturbations for both waves. 

At n(3' e ~ 1, the Vj3' e n term in Eq. ([76]) starts to play 
a role. The phase speed of the Alfven wave starts to in- 
crease in a way characteristic of the kinetic Alfven wave 
[32| and the character of the fluctuations (magnetic or 
acoustic) reverses for the two solutions. At very large re 
[for which our equations are probably not accurate (un- 
less /3i <C /3 e ) because they do not include standard FLR 
corrections, e.g., (k±pi) 2 = \ npi 1], the solution on 
the Alfven wave branch (solid curve) becomes an elec- 
tron sound wave with w/fcii = wth cJ the electron thermal 
speed. On the other hand, the solution on the sound 
wave branch has phase velocity equal to the sound speed 
based only on the ion temperature. 



VI. SIMULATION RESULTS 

We have implemented a two-dimensional reduced mag- 
netohydrodynamic (MHD) finite-difference simulation 
using the equations in this paper (see Table [TTJ) . One of 
the two dimensions of the simulation is the direction of 
the background magnetic field. The code uses generalized 
orthogonal coordinates [34[ , and is second-order accurate 
with respect to time and fourth-order accurate with re- 
spect to space. For this paper, we use an insulator bound- 
ary at the ends of the simulation encountered by mov- 
ing along the background magnetic field, and a hard- wall 
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FIG. 1: Plotted versus n (normalized k\ denned in Table U}, 
(a) the normalized squared phase speed V and (b) the ratio of 
relative magnetic fluctuations to relative density fluctuations 
SBy/Sn e = (SB y / Bo)/ (Sn e /no) for the linear dispersion re- 
lation Eq. (|76[l for a homogeneous two-fluid isotropic plasma 
with parameters defined in Table [I] 



TABLE II: Reduced fluid equations used in the 2-D reduced 
MHD simulation code. 



Continuity for each species 


Eq. 




Parallel momentum for each species 


Eq. 




Perpendicular pressure for each species 


Eq. 


§B 


Parallel pressure for each species 


Eq. 


(Ml) 


Reduced-fluid velocity for each species 


Eq. 




Reduced fluid displacement 


Eq. 


CD 


Effective potentials 


Eq. 


m 


Perturbed electric field 


Eq. 


gj 


Perturbed magnetic field 


Eq. 


© 


Reduced Poisson equation 


Eq. 




Reduced Ampere equation 


Eq. 


J531) 



perfect conductor boundary at the ends of the simulation 
encountered by moving within the simulation plane per- 
pendicular to the background magnetic field [35[ . These 
boundary conditions are energy conserving in the sense 
that there is no flux of energy out the boundaries. There 
is one modification of the equations that we made in the 
simulation code. We set b = 1 [defined in Eq. pTjl ] and 
dropped the magnetic term in pj_ [Eq. (jlip ] only within 
the magnetization current ()28|1 that appears in the re- 
duced parallel- Ampere equation ([5*1)1 . Making both of 
these changes together still maintains energy conserva- 
tion. The assumption is that the change in the total 
magnetic field amplitude caused by the (perpendicular) 
perturbation of the magnetic field is small, an assump- 



tion that is well satisfied for a low beta plasma such as 
occurs in the dipole magnetosphere at low altitudes. 

The main purpose of the simulations is to demonstrate 
good energy conservation, a major advantage of our La- 
grangian formulation. Note that the standard convective 
E x B and magnetic-flutter nonlinearities vanish in two 
dimensions (since ue • V = = Bj_-V), which leaves 
only the parallel convective term U||bo • V ^ 0. The non- 
linear FLR corrections, on the other hand, involve the 
differential operator 



X U E 



B_ 

Bo 



•V ? 0, 



which does not vanish in two-dimensional geometry. The 
use of a two-dimensional simulation geometry therefore 
enables us to focus our attention on the nonlinear FLR 
effects considered in the present work. 

The energy density £ is given in Eq. ([59]). We write 



with 



£ — £b + £i 

£b = 
£e = 
£k\\ = 

£k± = 



£ 



K\ 



■■K± 



+ £ 



£ 



£p|l = 



8^ 
2 ^ 

S 

2 ^ 

S 



m s n s 



P\\s 



Xs 



Ue + U|| 



P\\s0 



P±s0 



B_ 

'B 



(77) 

(78a) 
(78b) 
(78c) 

(78d) 

(78e) 

(78f) 



Note that the perturbation fields and are phys- 
ical fields, the CGL-pressure components (p± s ,P\\ s ) are 
physical fields in the standard zero-Larmor-radius limit 
[{p±sO,P\\so) denote initial (equilibrium) values], and the 
fields (n s ,uu s ) are gyrocenter-fluid fields related to their 
corresponding physical fields by Eq. (fT""f|) . This division 
of terms is somewhat different than that of Sec. IIII B[ 
where the energy density is divided into terms that best 
demonstrate the pathways of energy flow. Here, we chose 
to express the energy density terms in such a way that 
each term is positive definite [except for Eqs. (|78e[) - (|78fp ]. 
For the results shown in the following simulation plots, 
the energies are volume averaged over the system vol- 
ume [see Eq. ([To")) ]. Using the finite difference method, 
this means that we add up the products of the energy 
density and grid cell volume at each grid point. The re- 
sulting quantities are energies rather than energy density, 
but for the rest of this section, we use the same variable 
names defined in Eqs. ([75)1 (labels for the energy density 
terms) for the energy terms. 
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FIG. 2: Energy terms (as described in the text) versus time 
for a linear (5v v /Va = 10~ 10 ) simulation in straight geometry 
with parameters given in Table U We plot in (a) , all the 
energy terms except the change in the total energy 8£, in (b), 
5£ and the electric field (displacement current) energy £e, 
and in (c), just 8£. 



A. Linear simulation in straight geometry 

Figure [2] shows energy terms [based on energy densities 
defined in Eqs. (|78p] and the change in the total energy 
from the beginning of the run, 8£ , versus time t for a 
linear (Sv v /Va = 1CP 10 ) simulation in straight geometry. 
The parameters arc given in Table |TJ The energy terms 
are integrated over space and normalized so that the en- 
ergy of the constant equilibrium magnetic field Bq would 
be unity; time is normalized to an arbitrary normaliza- 
tion length (written in the figure caption as Re) divided 
by the Alfven speed. The code was initialized with a sinu- 
soidal (in both simulation directions) wave perturbation 
of the out-of-plane velocity u^; consistent with a wave on 
the Alfven branch in Fig. [1] (solid curve). For these pa- 
rameters, the total time of the simulation run was equal 
to the wave period of the wave determined from Eq. (|76|) 
on the Alfven branch [plus sign of quadratic equation for 
V in Eq. ([76])]. 

Two major results can be seen from Fig. [5] First, the 




4 Z(M 8 



FIG. 3: Grid points (dots) of the dipole simulation plotted 
in real (Cartesian) coordinates X and Z (direction of dipole 
axis). [Only every fourth point is plotted in the Q (parallel) 
direction and every sixteenth point in the L (across field) di- 
rection.] The curvilinear coordinate directions are also shown, 
with the parallel coordinate Q varying along the magnetic 
field and the coordinate L varying across the magnetic field. 
Also shown are contours for the initial perturbation in $ cen- 
tered at (X,Z) = (7,0). 



resulting oscillations are consistent with the wave pe- 
riod from Eq. (|76p since the run has two complete os- 
cillations of the (quadratic) energy terms. Secondly, the 
change in the total energy (labeled 6£) is much smaller 
than the change in any individual energy term as can 
be seen by comparing the size of the terms in Fig. [5^, 
b, and c. Whereas rough constancy of the total energy 
(Fig. [2^) does not well demonstrate energy conservation, 
the fact that the change in the total energy is smaller 
than the change in energy of the individual terms (Fig. [2^ 
and b) does usually indicate good energy conservation. 
It shows that the error in the energy is smaller than the 
energy associated with the dynamics of the particular 
term. Fig. [2t> shows that the change in the total energy 
is even less than the energy of the electric field associated 
with the displacement current (usually neglected). In 
Sec. IVICI we present a more detailed convergence study 
demonstrating the quality of the energy conservation. 



B. Linear simulation in dipole geometry 

Next, we show linear results in dipole geometry. The 
simulation grid and coordinate system are described in 
Fig. [3J We use a 256 x 256 grid in the coordinates Q and 
L (Fig. [3]). The coordinate L is the L shell used in mag- 
netospheric physics, and is equal to the radial distance 
at the magnetic equator (Z = in Fig. [3]) in units of the 
Earth radius Re- The system is perturbed with a peaked 
distribution of $ at L = 7 at the magnetic equator (con- 
tours in Fig. [3]). The parameters of this two-fluid plasma 
simulation are given in Table [J (where k is defined at the 
magnetic equator based on the scale length of the initial 
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value for $). 

Figure [5] shows contours of El, An, SB V (out of 
a plane component), and Sn e with respect to the curvi- 
linear coordinates Q and L. These arc plotted at four 
different times normalized to Re /Vao, where Vao is the 
Alfven speed at L = 7 and at the magnetic equator. The 
initial perturbation in <f> is plotted in the top left panel 
(t = 0) . The shape of the contours is different from those 
of Fig.[3]because of the different coordinate system, but it 
is clear from both figures that we have a single monotonic 
peak. In addition to the initial perturbation in <£>, there 
is also an initial perturbation in n e (bottom left panel). 
Not shown are the initial perturbations for n,;, p\u, pj_%, 
pn e , and p± e . (There is no initial magnetic field pertur- 
bation.) All the variables are initialized with parameters 
consistent with a linear wave with parallel and perpen- 
dicular wavelengths corresponding to the scale lengths of 
the initial perturbation in the parallel and perpendicular 
directions. 

As can be seen from Fig. [H the initial perturbation 
breaks up into two traveling waves that move along the 
magnetic field (left or right). The traveling waves have 
a magnetic perturbation. The Alfven speed ~ Bq and 
Bo ~ L~ 3 . Because of this, the perturbation travels 
faster at low values of L. In addition, the real length 
of the field line is less at lower L (Fig. [3]). Therefore 
the perturbation travels much faster along the curvilin- 
ear coordinate Q at the low values of L. There are two 
lobes of the perturbation El or 5B V . By the final time 
of the simulation, the lobes at lower L ~ 6.5 have al- 
ready reflected off of the insulating boundary at Q = ±1, 
while the lobes at higher L ~ 7.5 have not. In the MHD 
limit, one would associate the density perturbation with 
a sound wave rather than an Alfven wave, but Fig. 0] 
shows that all the perturbations travel together, consis- 
tent with a coherent linear wave. 



C. Nonlinear simulation in dipole geometry 

Finally, we ran a third simulation with the same pa- 
rameters as for the linear run in dipole geometry (Fig.O, 



but with initial velocity perturbation Sv v /Va = 0.3. This 
amplitude is quite large, and when a simulation is run in 
straight geometry with a purely sinusoidal perturbation, 
the fluctuations of the energy terms are not regular as 
they are in Fig. [2] (not shown; especially irregular are the 
fluctuations in the parallel kinetic and pressure energies). 
Figure [5] (similar to Fig. [21) displays the energy terms for 
the nonlinear dipole simulation, and demonstrates that 
the total simulation energy is well conserved in this case 
also. 

To demonstrate the convergence properties of this non- 
linear energy conservation, we now focus our attention on 
the simulation time interval ranging from t (Vao/ Re) = 
to 1 in Fig. [5l where the largest changes in energy are 
for the magnetic energy and the perpendicular kinetic 
energy. Figure [6] shows the log-log plot of the error in 
the total energy 5£ normalized to the magnetic energy 
£b (which is approximately the same for all the simula- 
tion results) as a function of the normalized time step 
At (Vao/ Re) for different grid resolutions N t = 128, 256, 
and 512 (the same in both the field-line Q-direction and 
radial L-direction). In the best-case scenario, we should 
see the following behavior as we reduce At at a fixed 
value of N{: Since the time step algorithm is a second- 
order predictor-corrector scheme (leapfrog trapezoidal), 
the error in the energy should go down as a factor of four 
for every reduction in At of a factor of two. The error 
in the energy should then converge to a constant value 
limited by the spatial resolution (as observed in Fig. [(5J) ; 
note that if we kept decreasing At, the error in the energy 
would eventually rise because the computer calculations 
would not have the precision necessary to accurately solve 
the equations. Because the spatial-differencing scheme 
for our simulation code is spatially fourth-order accurate, 
we should see a decrease in the energy error of a factor of 
2 4 = 16 each time we double the resolution (for a fixed 
time step At), and we do in fact see that decrease in the 
time-resolved error in Fig. [5] as we increase the number 
of grid points N{ from 128 to 256 (decrease in error of 
2.1 x l(T 6 /i-3 x 1(T 7 ~ 16) and from N, = 256 to 512 
(decrease in error of 1.3 x 10~ 7 /8.0 x 10~ 9 ~ 16). 



While we regard a detailed description of the physics 
of the new FLR terms as beyond the scope of this paper, 
we can easily demonstrate that they have an appreciable 
effect on the energy conservation. If we run a simula- 
tion with Ni = 256 and At(V A o/RE) = 3.125 x 1CT 4 
(converged with respect to time for this grid resolu- 
tion), but setting K p in Eq. (|68|) equal to zero (that 
is, setting the nonlinear ponderomotive force density 
equal to zero), we find the normalized energy error to 
be 5£/£b = 1-4 x 10 -3 , which is much larger than the 
value 1.3 x 10~ 7 shown in Fig. [51 The fact that this 
error is less than unity indicates that the zeroth order 



physics of Alfven waves (energy transfer between mag- 
netic and perpendicular kinetic energy) is being correctly 
described. However, the ratio of the change of the total 
energy is 43 times greater than that of the electric field 
energy and 1200 times greater than the parallel kinetic 
energy, showing that the parallel dynamics and dynamics 
associated with the displacement current are not at all 
well described. 



13 



t = OR/V AQ t=1R/V AQ t = 2R/V M t = 3R/V A0 



L 




6 r,,,,i,,,,T,,,,i,,,,T,,,,i,,,,T,i,,i,,,i" 
-1 1-1 1-1 1-1 1 

Q 



FIG. 4: Contours of (from top to bottom) <E>, El, An, SB^,, and Sn e at four different times (arranged from left to right) 
indicated at the top of the plot. In each panel, the contours are plotted versus the curvilinear coordinates Q and L and positive 
(negative) contour levels are solid (dashed). 



VII. SUMMARY 



The variational derivation of the reduced fluid equa- 
tions has revealed the existence of a new type on nonlin- 
earity in reduced fluid dynamics. While standard fluid 
nonlinearities appear in the convective derivative oper- 
ator u-V (e.g., E x B and magnetic- flutter nonlinear- 
ities), the new nonlinear terms presented here can be 
described as nonlinear FLR effects that appear as cor- 
rections (ITBl of the electromagnetic potentials $ and A\\. 
These nonlinear FLR-corrected potentials are then used 
to construct the magnetic and electric fields ([?!?)) and 
(|75p that appear in the parallel reduced-fluid equation 
(I74p . which also contains a ponderomotive- force correc- 



tion (frj§)) to the standard CGL-pressure force density. 

The linear properties of the equations for a two-fluid 
homogeneous plasma were described, and linear and non- 
linear simulations demonstrated that the equations de- 
scribe the coupled dynamics of Alfven and sound waves 
and that the simulation energy is conserved in both 
straight and dipole geometry. 

Lastly, we note that the single limitation on the back- 
ground magnetic field in the present work was associ- 
ated with the absence of parallel current along the field 
lines (i.e., bo • V X bo = 0). This constraint was used as 
a simplifying assumption [see Eqs. (|46|) - (|48)) ] associated 
with the perturbed magnetic field V X (A\\ bo) = 
having no component along the background field lines 
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FIG. 5: Energy terms versus time for a nonlinear (Sv^/Vao = 
0.3) simulation in dipole geometry with parameters given in 
Table HI We plot in (a), the largest energy terms, in (b), the 
smaller energy terms with the change in the total energy 8£, 
and in (c), just 5£. 



APPENDIX A: PUSH-FORWARD 
REPRESENTATION OF FLUID MOMENTS 

The purpose of this Appendix is to present the the- 
oretical foundations [2(| that establish the relation be- 
tween fluid moments of a particle distribution function 
-Fphy S (z,t) in particle phase space and fluid moments of 
a reduced distribution function F(Z, t) in reduced phase 
space. Here, the near-identity (reversible) transforma- 
tion T t : z — * Z from particle phase space to reduced 
phase space is represented in asymptotic form as 



1 



where e <C 1 is a small ordering parameter and 
(Gi,G2,---) represent the Lie-transform vector fields 
that generate the transformation T t . Furthermore, this 
transformation induces the pull-back operator T £ : F — > 
Fphys — T e -F and push-forward operator T7 1 : F p \ lys — > 
F = T^Fphys between phase-space functions, which 
both preserve the scalar invariance property F p ^ ys (z, t) = 
F(Z,t). 

The reduced displacement 



t; 



(Al) 



(e.g., gyroradius) between the reduced position (e.g., 
guiding-center position) and the push-forward of the par- 
ticle position is of particular importance here and it is 
expressed as 

1 



-eGj - e G 



Gi • dGn + 



in terms of the Lie-transform generating vector fields 
(Gi,G 2 ,---). 



General reduced fluid moments 



(bo-B^ = 0). More general magnetic geometries with 
b • V X b 7^ (e.g., tokamak geometry) can also be 
treated within a variational formulation and will be the 
subject of future work. 
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We begin by considering an arbitrary fluid moment on 
physical (phys) particle phase space: 



("M)pbyB = J d 3 pxF phys 



= J d 6 z X <5 3 (x-r)F phys 

= Jd 6 ZT- 1 X 5 3 (X + p £ -r)F,(A2) 

where x is an arbitrary function on particle phase space 
and T~ 1 x is its push-forward on reduced phase space. 
Here, [x]ph ys denotes the physical particle-momentum av- 
erage of X with respect to Fpiiys an d ^phys is the particle 
fluid density. Upon Taylor expanding Eq. (|A2D in powers 
of p € and integrating by parts, we find the nonlinear FLR 
expansion 



(n[ X ]) phys = n [T^x] - V- (n [p^x]) 



(A3) 
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FIG. 6: For the nonlinear simulation in dipole geometry (for the time interval t(VAo/RE) — to 1 in Fig. [5}, the logarithm 
of the normalized total energy 8£/£b is shown as a function of the logarithm of the normalized time step At (Vao/Re) for 
different numbers of grid points Ni = 128, 256, 512 (the same number in both directions). For a fixed number of grid points 
Ni (lines are used as guides), each successive normalized time step is reduced by half as we proceed to the left. In addition, 
the impact of setting K p in Eq. (|68[) equal to zero (that is, setting the nonlinear ponderomotive force density equal to zero) is 
reflected by a large jump (four orders of magnitude) in the total-energy non-conservation. 



where [• • • ] denotes the reduced-momentum average with 
respect to F and n denotes the reduced fluid density. 
Note that the push-forward representation (|A3|) is re- 
versible and its inverse yields the pull-back representa- 
tion [5j 

n[x] = "phys [TeXlphys + V< ("phys [Pe T eX] p hys ) 

+ ■ • • , (A4) 

where \ is an arbitrary function on reduced phase space 
and T e x is its pull-back on particle phase space. 

We now consider the fluid moments associated with 
fluid density (x = 1) and fluid velocity (\ = v = dx/dt). 
First, the push-forward representation for the particle 
fluid density 



^phys 



[Pe] ~ -V-(nM)| • (A:,) 



where we have retained the second-order ( "quadrupole" ) 
term, which will prove useful in what follows [201 ] . 

Next, using the definition (|A1|) . we consider the push- 
forward representation for the particle fluid velocity 



t; 



e dt e 



(T £ - 



~dT 



d e p e 
dt 



(A6) 



where d e /dt denotes the reduced Vlasov operator so 
that deX/dt denotes the reduced particle velocity (e.g., 
guiding-center velocity) and d e p e /dt denotes the re- 
duced "displacement" velocity (which may include the 



gyroangle-dependent perpendicular velocity and the 
gyroangle- independent polarization velocity). The push- 
forward representation for the particle flux Tp^ = 
in [v]) phys : 



■ phys 



d e X 
~dT 

+ V X 

r + r pol - 



d_ 



Pe X 



1 d e p e 



2 dt 



+ 



d e X 
~~dT 



(A7) 



is expressed in terms of the reduced flux T = n [d £ X/dt], 
the polarization flux 



- pol 



d_ 

di 



n [p £ ] , 



(A8) 



and the reduced (divergenceless) magnetization flux 



V X 



Pe X 



deP e 

dt 



Pe X 



cLX 



dt 



which is itself decomposed in terms of an intrinsic contri- 
bution (first term) and a moving-dipole contribution (sec- 
ond term) [25]. Note that the correct derivation of the 
intrinsic contribution relied on keeping the quadrupole 
contribution in Eq. (| A5|) . 

The push-forward representations (| A5|) and (|A7|) pre- 
serve the conservation law of particles through the con- 



16 



tinuity equation 

t^phys 







at 



d_ 

at 



+ V • Tphys 

V- (n[pj) 



v • (r + r 



pol 



■ mag; 



dn 

where the reduced polarization effects cancel each other 
exactly and the reduced magnetization flux is divergence- 
less (by definition). 

Lastly, using these push-forward representations (|A5|) 
and (|A7|) , we may write the push- forward representation 
for charge density as 

Pphys = S^phys = P - V-P e , 

where p = ^ qn and the reduced polarization vector is 



~P e = qn [p e ] 



and the push-forward representation for current density 
as 

J phys = <7r phys = J|-^ + cVxM e , 
where J = q T and the reduced magnetization vector 



qn 
c 



Pe X 



1 d e p e 

2 dt 



~dT 



The reduced Maxwell equations are, thus, expressed in 
terms of the "macroscopic" fields D = E + 47r P e and 
H = B - 4?rM f as 



V • D = 4tt p and V X H - 



1 3D 

c ~dt 



47T 



2. Guiding-center and gyrocenter fluid moments 

Let us now obtain explicit expressions for the reduced 
fluid moments associated with the guiding-center and gy- 
rocenter phase-space transformations. 

We begin with the guiding-center (gc) transformation 
for which p t = p (and we ignore electric and magnetic 
perturbation fields until we consider the next transfor- 
mation). By substituting the gyroangle-dependent gyro- 
radius p in Eq. (|A5|) . we easily find that the physical 
and guiding-center densities are equal rt p h ys = n gc , since 
we ignore standard FLR corrections in the present work; 
in the same spirit, the perpendicular and parallel pres- 
sures are also identical in the physical and guiding-center 
fluid descriptions. The relation (|A7|) between the physi- 
cal particle flux and the guiding-center particle flux, on 
the other hand, yields the expression 



- phys 



r gc - vx 



Pa 



(A10) 



where the guiding-center flux is 



r KC = nuu b + — — X Vlnl? + p\\ b- Vb 

11 mil \ 

only the guiding-center parallel magnetization 
(— p±b/mfl) survives gyroangle averaging, and V-P 
denotes the divergence of the CGL pressure tensor (Til)j) . 
Here, the parallel fluid velocities are also identical in the 
physical and guiding-center fluid descriptions (labeled 
tin). Note that this expression is equivalent to the 
reduced-fluid velocity (|4"T)) in the absence of electric and 
magnetic perturbation fields. 

Next, we consider the gyrocenter (gy) transformation 
(now involving the electric and magnetic perturbation 
fields) for which we take [p e ] = p ± . In this case, using 
Eq. (|A5|) . the guiding-center (and physical) density is 
expressed 



V • (n gy p A 



(All) 



The relation (|A7[) between the physical particle flux and 
the guiding-center particle flux, on the other hand, yields 
the expression 



d_ 

dt 



( n gy Pi 



V X 



(n«||)gy p ± X b 



(A12) 



which includes the polarization-drift flux and the 
(moving-dipole) magnetization flux to first order in 
NFLR corrections. This expression can be used to de- 
rive the parallel flux relation 



llgc 



(nun) 



bn • V X 



(nu\\)gy p ± x b 



"gy 



V> (%y P±)] u ||gy 

(Pj_-VV||gy) + 



(A13) 



where we have ignored the effects of background magnetic 
curvature in obtaining the second equality. By using the 
relation (jAllj) between the guiding-center and gyrocenter 
densities, we now obtain the relation 



llgv 



• Viti 



gy 



(A14) 



where we have ignored higher-order NFLR corrections. 

Lastly, we combine the guiding-center and gyrocenter 
push-forward relations to obtain the expressions for the 
physical (particle) density 



"-phys = n - V- (n px) , 
and the physical (particle) flux 

Iphys ^ ^ 

+ V X (np x XU| b 



(A15) 



Wt (n P± ) 



(A16) 



in terms of the gyrocenter density n and the gyrocenter 
velocity u (used in the text). 
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APPENDIX B: REDUCED-FLUID 
POLARIZATION DISPLACEMENT 



where we neglected curvature effects in going from the 
first line to the second line. Identifying 



Since p^ is the key quantity describing our new inter- 
pretation of the reduced fluid equations (first presented 
in Ref. [3), it is important to have a good understanding 
of what it is. The displacement vector p^ is shown here 
to possess a s imp le interpretation as the time integrated 
inertial drift [361 ] . 

Neglecting curvature, we can easily derive a formula for 
p^ from the single particle equation of motion. Writing 
the particle velocity as 



V + Vij_, 



v = u\\ b 



Bi 



Ue, 



we approximate the equation of motion as 



^ v o Q I ,, 1 

— - ~ — I E + -v X B 

at m 

= — Vi ± X B . 



Solving for Vi_i_, we find 



Vi_L 



b dv 

Uq at 

d ( 1 u 

dt u bo x vo 



d_ / b^ 
dt 1 n 



X U E 




(Bla) 
(Bib) 



(B2) 



(B3) 



Vl. 



dp ± /dt 



(B4) 



immediately leads to Eq. (fTT|) . 

Next, we now present a simple derivation of the NFLR- 
corrected potentials (fT3|) based on the electromagnetic 
one-form A = A • dx — <f> cdt. By expanding the NFLR- 
corrected one-form 

A p = A(x + p ± ,t)-(dx + dp ± ) - $(x + p ± ,t) cdt 

to first order in p ± , we obtain 

A p = A + p ± - (VA-dx - V$ cdt) + A ■ dp ± 
= (A - p ± x B) • dx - ($ - p ± ■ E) cdt 
+ d(A-p x ) 

= (A, + Vt?) - dx - Up - ~|H cdt, (B5) 



where 



/ $ - p ± • E 



(B6) 

\A - p ± xB 

and r\ = A • p x is treated as a gauge term. In Eq. (|13[) . 
the vector potential is A = Ay b (i.e., ry = 0) and A\\ p 
b • A p . 
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